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Abstract 

We  discuss  issues  related  to  modeling  of  nonlinearities  and  hysteresis  arising  in  a  class  of 
magnetorheological-based  smart  elastomers.  The  dynamic  models,  intended  for  use  in  parameter 
estimation  and  control  problems,  are  presented  in  the  context  of  simple  elongation  of  a  filled 
rubber-like  rod.  Theoretical,  computational  and  experimental  results  are  given. 


1  Introduction 

Smart,  material  structures  and  fluids  are  generally  understood  to  be  structure  and  fluid  composites 
that  possess  the  capability  to  sense  and  actuate  in  a  controlled  manner  in  response  to  variable  ambi¬ 
ent  stimuli.  These  are  in  actuality  smart  material  systems  which  involve  combinations  of  advanced 
sensors,  actuators  and  microprocessors.  Effective  practical  use  of  these  systems  in  specific  appli¬ 
cations  depends  on  fundamental  developments  related  to  a  number  of  important  modeling  issues 
involving  these  composites.  In  particular  one  must  have  (i)  models  for  the  composite  host  system 
including  sensors  and  actuators  and  (ii)  models  describing  host  system  responses  to  input  signals 
to  the  actuator.  The  first  of  these  is  the  inactive  host  system  model  while  the  latter  refers  to  the 
system  undergoing  actuation  in  response  to  stimulation.  For  example,  detailed  discussions  of  issues 
related  to  models  for  self-sensing,  self- actuating  structures  based  on  piezoceramic  sensors/actuators 
can  be  found  in  [9].  Piezoceramic-based  smart  systems  are  a  rapidly  maturing  technological  field 
with  a  large  literature  (a  substantial  number  of  references  are  given  in  [9]).  Another  class  of  systems 
based  on  magnetorheological  (MR)  solids  and  fluids  is  far  less  developed  and  questions  related  to 
these  systems  are  the  main  focus  of  our  discussions  in  this  presentation. 

MR  elastomers  [26], [27],  which  are  solid  analogs  of  MR  fluids,  are  rubber-like  composite  struc¬ 
tures  filled  with  active  as  well  as  inactive  substances.  Magnetically  permeable  particles  (such  as 
carbonyl  iron)  are  added  to  a  viscoelastic  polymeric  material  prior  to  crosslinking.  A  strong  external 
magnetic  field  is  applied  before  and  during  crosslinking.  This  field  induces  dipole  moments  within 
the  particles,  which  seek  minimum  energy  states.  Particle  chains  with  collinear  dipole  moments  are 
formed  and  curing  of  the  polymeric  host  material  locks  the  chains  in  place.  The  resulting  product 

‘Invited  Keynote  Lecture,  4th  International  Conference  on  Optimization:  Techniques  and  Applications,  Perth, 
Australia,  July  1-3,  1998. 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

JUL  1998 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-1998  to  00-00-1998 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Estimation  and  Control  Related  Issues  in  Smart  Material  Structures  and 

5b.  GRANT  NUMBER 

r  iuius 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

North  Carolina  State  University, Center  for  Research  in  Scientific 
Computation, Raleigh, NC, 27695-8205 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

17 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


is  a  composite  material  with  variable  elastic  modulus  (maximum  modulus  changes  of  between  30% 
and  40%  have  been  reported  in  experimental  data  [26]  in  response  to  an  applied  external  magnetic 
field). 

Models  of  the  response  of  the  elastomer  to  an  applied  magnetic  field  involve  induced  stress  as  a 
function  of  induced  magnetic  flux  density.  This  in  turn  involves  understanding  the  dependence  of 
the  elastic  modulus  on  the  induced  magnetic  flux.  Modeling,  which  is  truly  in  its  infancy,  can  be 
based  on  magnetic  dipole  interactions  between  adjacent  particles.  Since  the  value  of  the  magnetic 
permeability  varies  dramatically  between  particles  and  host  material,  it  is  to  be  expected  that 
homogenization  techniques  [10],  [16],  [22]  will  play  a  significant  role  in  any  careful  modeling  of 
effective  moduli.  Early  results  [32]  suggest  that  the  magnetic  permeability  is  a  nonlinear  function 
of  the  magnetic  potential  and  hence  nonlinear  homogenization  [15]  formulations  will  be  necessary 
to  develop  models  in  category  (ii)  above. 

Our  discussions  here  will  be  on  models  for  the  dynamical  response  of  the  composite  host,  i.e., 
inactive  filled  elastomers  or  filled  viscoelastic  structures.  Thus  we  focus  on  issues  related  to  require¬ 
ment  (i)  above.  Such  models  involve  nonlinear  and  hysteretic  formulations  in  a  significant  way.  We 
describe  some  of  our  efforts  with  model  development,  estimation  and  experimental  verification. 


2  Nonlinear  Constitutive  Laws 


Our  early  efforts  focused  on  understanding  the  nonlinearities  inherent  in  the  dynamic  response 
of  rubber-like  compounds  to  applied  loads  and  impulsive  disturbances.  First  we  outline  the 
development  of  our  basic  model,  then  turn  to  theoretical  foundations  and  experimental  results. 

Most  models  for  elastomers  found  in  the  literature  are  based  on  strain  energy  function  (SEF) 
and  finite  strain  theories  (see  [6]  and  references  there).  To  illustrate  our  approach  (which  is  based 
on  the  same  principles),  we  take  a  simple  example  (for  a  detailed  development  see  [6]):  an  isotropic, 
incompressible,  rubber-like  rod  with  a  tip  mass  undergoing  simple  elongation  with  a  finite  applied 
stress  in  the  principle  axis  direction  aq  =  x.  For  neo-Hookean  materials  the  finite  stress  theory 
(or  the  Mooney  theory  with  SEF  U  =  C\{I\  —  3))  leads  to  a  true  stress  T  =  y(Af  —  ^j-),  or  an 
engineering  stress 


S-L-&X 

Here  A*  are  the  principle  extension  ratios  which  represent  the  deformed  length  of  unit  vectors 
parallel  to  the  principal  axes  (the  axes  of  zero  shear  stress),  and  the  first  strain  invariant  I\  is 
I\  =  A^  +  A|  +  A|.  In  this  case  the  dynamical  problem  reduces  to  a  1  —  D  problem.  In  terms  of  the 
finite  strain  err  and  the  deformation  u  in  the  x  direction  we  have 
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This  can  be  used  in  the  Timoshenko  theory  for  longitudinal  vibrations  of  a  rubber  bar  with  a  tip 
mass  to  obtain  (p  =  mass  density,  F(t)  =  applied  external  force,  Ac  is  the  cross  sectional  area,  M 
is  the  tip  mass,  g  is  the  gravitational  constant) 
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where  S,  the  internal  stress  resultant,  is  given  by 
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Here  E  is  the  generalized  modulus  of  elasticity  and  </(£)  =  1  +£  —  ( 1  +  £  )-2.  A  Kelvin- Voigt  damping 
term  is  included  in  the  stress  as  a  first  attempt  to  model  damping  ( G'/j  is  the  Kelvin- Voigt  damping 
coefficient).  This  leads  to  the  nonlinear  partial  differential  equation 
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v=£  +  E(t )  +  Mg  (4) 


«(/,  0)  =  0,  n(0,  x)  =  A(x),,  •«,(().  x  )  -  0.  (5) 

for  dynamic  longitudinal  displacements  of  a  neo-Hookean  material  rod  in  extension.  Here  the  initial 
configuration  is  given  by  A(i)  and  g(()  =  £  —  (1  +  £)-2.  This  model  in  variational  formulation  is 
given  by 

pAciitt  +  Aiu  +  A2ut  +  D*g(Du)  =  F  in  V* ,  (6) 

where  V  is  an  appropriately  chosen  Hilbert  space. 

Several  theoretical  and  practical  questions  arise  in  connection  with  this  model.  The  most 
fundamental  is  the  well-posedness  of  the  initial  boundary  value  problem  associated  with  (3)-(5).  ft 
is  shown  in  [4]  that  a  unique  weak  solution  of  ( 3 ) - ( 5 )  exists  even  under  more  general  circumstances, 
i.e.,  for  a  broad  class  of  nonlinearities  g.  This  is  important,  since  comparison  between  experimental 
data  and  preliminary  numerical  calculations  suggested  that  the  neo-Hookean  nonlinearity  is  not 
adequate  to  describe  the  behavior  of  filled  elastomers.  To  estimate  the  correct  g  and  the  other 
unknown  parameters  in  (3)  (5),  namely,  p,  K.  and  C#,  free  release  experiments  were  performed 
(for  details  see  [6]). 

For  the  experiments  a  slender  rod  composed  of  unfilled  natural  rubber  was  used.  The  rod  had  a 
tip  mass  to  guarantee  that  it  remains  in  simple  extension  and  compression  is  not  present.  Initially, 
the  rod  was  lifted  so  that  it  was  at  its  natural  length.  Then  the  support  was  removed,  allowing 
the  mass  to  fall  freely.  The  rod  achieved  approximately  34%  maximum  dynamic  strain  during  this 
test.  The  data  collected  by  the  load  cell  on  top  of  the  sample  were  used  to  estimate  the  unknown 
parameters. 

In  one  general  parameter  estimation  formulation,  equation  (6)  takes  the  form 

pAciitt  +  A\(q)u  +  A2{q)ik  +  D*g(q)(Du )  =  F(q)  (7) 

•*(())  -  A(;r),  «;(0)  -0.  (8) 

where  the  structural  operators  A\,A2,  the  nonlinearity  g  and  the  input  F  have  all  been  parame¬ 
terized  by  a  vector  (possibly  infinite  dimensional)  parameter  q  that  must  be  estimated.  Here  the 
parameter  q  takes  values  from  an  admissible  parameter  set  Q.  Suppose  that  we  have  a  set  of 
measured  observations  z  =  ,  corresponding  to  measurements  (e.g.,  displacements,  velocities) 

taken  at  time  tj.  In  a  general  least  squares  parameter  estimation  problem,  we  seek  to  minimize  the 
least  squares  output  functional 


J{q,z) 


C2  {C\{ii(ti,-,q)}-  { Zi }}  2 


(9) 


over  g  €  Q,  where  {u(  /,,  •;  q)}  are  the  parameter- dependent  solutions  of  (7  )  (8)  evaluated  at  time 
ti,  i  =  1,2,..  .K ,  and  |  •  ]  is  an  appropriately  chosen  Euclidean  norm.  The  operators  G'i,  C'2  depend 
on  the  type  of  the  collected  data.  For  example,  if  Z{  is  time  domain  displacement,  velocity  or 
acceleration  at  a  point  x ,  then  C\  involves  differentiation  (0,1  or  2  times,  respectively)  with  respect 
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to  time  followed  by  pointwise  evaluation  in  t  and  x.  The  operator  C2  is  the  identity  in  the  case 
of  time  domain  identification,  while  it  is  related  to  the  Fourier  transform  if  we  consider  fitting  the 
data  in  the  frequency  domain  (see  Chapter  5  of  [9]  for  details). 

In  this  formulation  the  minimization  problem  involves  an  infinite  dimensional  state  space  and 
(in  general)  an  infinite  dimensional  admissible  parameter  set  Q.  To  overcome  this  difficulty  and 
to  obtain  a  computationally  tractable  method,  we  use  the  general  ideas  described  in  [9].  Namely, 
let  Hn  be  finite  dimensional  subspaces  of  the  state  space  H ,  and  QM  be  a  sequence  of  finite 
dimensional  sets  approximating  the  parameter  set  Q.  Denote  the  orthogonal  projections  of  H 
onto  Hn  by  PN .  One  can  formulate  a  family  of  approximating  estimation  problems  with  finite 
dimensional  state  spaces  and  finite  dimensional  parameter  sets  in  the  following  way:  find  q  G  QM 
which  minimizes 
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where  uN(t. ;  q)  G  HN  is  the  solution  to  the  finite  dimensional  approximation  of  (7)-(8)  given  by: 


v,  4>)v*,v  +  1$ ,  4>)v*,v2  +  ( g(q)(DuN ),  Dcj)) 

=  {F{t\q)A)v*y2  (11.) 

uN( 0)  =  PN A(ar),  «f( 0)  =  0  (12) 


for  all  (f>  G  H N .  Here  V  c— ►  V2  ^  H  in  a  Gelfand  quintuple  (see  [9]). 

Solution  of  these  approximate  estimation  problems  ( 10 )-( 12 )  provides  one  with  a  sequence 
of  parameter  estimates  {qN,M}.  A  crucial  question  is  if  one  can  guarantee  that  this  sequence  (or 
some  subsequence)  converges  to  a  solution  of  the  original  infinite  dimensional  parameter  estimation 
problem.  Under  certain  suitable  assumptions  on  the  approximating  spaces  H  N  and  approximating 
sets  QM,  it  is  shown  in  [11]  that  qN’M  — ?  q*:  where  q*  solves  the  parameter  estimation  problem  for 
(9). 


The  above  scheme  was  implemented  using  linear  splines  to  generate  the  approximating  state 
spaces  and  to  approximate  the  nonlinear  term  g  (see  [6]).  Figure  1  depicts  the  best  fit  to  the  data 
that  can  be  achieved  using  a  linear  g,  and  definitely  suggests  that  a  nonlinear  function  is  required 
to  successfully  model  the  behavior  of  the  rubber  rod.  Figure  2  depicts  the  fit  achieved  using  a 
four-term  piecewise  linear  g.  This  fit  is  much  better  and  captures  the  nonlinear  behavior  as  well. 
However,  as  the  same  experiments  and  calculations  were  repeated  for  medium  and  highly  filled 
elastomer  rods,  the  best  attainable  fit  deteriorated,  indicating  that  in  these  cases  the  hysteretic 
behavior  also  must  be  taken  into  account. 


3  Hysteretic  Models 

3.1  Conceptual  Issues 

The  efforts  described  above  on  nonlinear  aspects  of  rubber  dynamics  also  conclusively  demonstrated 
the  (not  surprising)  importance  of  hysteresis  in  our  understanding  of  filled  host  dynamic  response. 
This  is  an  often  studied  and  seldom  resolved  aspect  of  viscoelastic  material  dynamics  which  usually 
begins  with  constitutive  stress-strain  formulations.  A  huge  literature  on  modeling  of  viscoelastics 
and  rubbers  currently  exists  (e.g.,  [12],  [14],  [19],  [34],  among  the  many  discussed  in  [35]).  Two 
types  of  models  for  stress-strain  relationships  can  be  found  in  the  literature  on  viscoelastic  mate¬ 
rials.  One  is  developed  on  the  basis  of  the  phenomenological  mechanical  behavior  of  the  samples, 
while  the  other  concentrates  on  the  microscopic  behavior  of  the  fibers,  such  as  the  changes  in  the 
crosslinking  and  the  contour  length  of  fibers,  and  on  the  different  relaxation  times  of  fibers.  The 
most  fundamental  model  is  the  Boltzmann  integral  model,  which  attempts  to  capture  the  viscosity 
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Linear  fit  to  load  cell  cdata  for  unfilled  rubber  sample,  small  weiaht  Linear  fit  to  load  cell  cdata  for  unfilled  rubber  sample,  small  weiaht 


Time  (seconds)  Frequency  (hz) 

Figure  1:  (left)  Time  domain  approximation  with  a  linear  g ,  and  (right)  the  FFT  of  the  solution 
and  the  data. 


Four  term  fit  to  load  cell  data  for  unfilled  rubber  sample,  small  weiaht  Four  term  fit  to  load  cell  cdata  for  unfilled  rubber  sample,  small  weiaht 


Time  (seconds)  Frequency  (hz) 

Figure  2:  (left)  Time  domain  approximation  with  a  four-term  piecewise  linear  g,  and  (right)  the 
FFT  of  the  solution  and  the  data. 
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of  the  material  and  the  history  dependence  of  the  stress  on  the  strain  and/or  strain  rate;  the  latter 
characteristic  is  a  signature  of  rubber  materials.  As  is  well  known,  the  Boltzmann  integral  can  be 
reduced  easily  to  some  well-known  differential  models,  e.g.,  Kelvin- Voigt  and  Maxwell.  Due  to  the 
fundamental  importance  of  the  Boltzmann  integral,  we  digress  here  to  briefly  outline  its  foundation. 

The  major  assumption  (a  principle  of  superposition)  made  by  Boltzmann  is  that  linear  models 
can  be  added  together  to  generate  more  complicated  models  (see  [28]).  This  assumption  is  referred 
to  as  the  Boltzmann  Superposition  Principle  (see  [18],  p.  6). 

One  can  generalize  the  basic  Kelvin- Voigt  and  Maxwell  models  (see  [35],  Chap.  3  for  a  detailed 
discussion)  to  relate  the  stress  and  strain  by 

s(t)  =  fT0./(/)[// ( /)]  where  a  =  (j0[H{t.)\  (13) 

a(t)  =  s0 Y(t)[H(t)]  where  s  --  c0[// ( / )]  (14) 

for  functions  J(t)  and  Y(t)  that  are  termed  the  creep  compliance  function  and  the  relaxation 
modulus  function,  respectively.  Here  H  is  the  usual  Heaviside  function.  Boltzmann  generalized  the 
above  models  to  account  for  variables  a(t)  and  e(t )  (e.g.,  the  results  from  stretch  and  relaxation 
tests)  by  considering  a  succession  of  infinitesimal  steps  ds(t )  for  (13)  or  a  succession  of  infinitesimal 
steps  da(t)  for  (14). 

As  a  result,  one  obtains 

da(t)  =  ds(s)Y(t  —  s)[H{t  —  s)]  .  (15) 

One  can  rewrite  ds(s)  =  dej^  ds ,  set  t,  >  s,  and  integrate  (15)  from  —  oo  to  t  to  obtain 

a (t)  =  I  Y(l  -  s)^j^-ds  . 

A  similar  expression  for  e  in  terms  of  the  kernel  ./  and  the  history  of  the  stress  rate  can  be  obtained. 
The  above  integral  is  referred  to  as  a  Boltzmann  integral  since  it  conforms  with  the  fundamentals  of 
superposition  as  enunciated  by  Boltzmann.  Physical  models  that  consist  of  a  continuum  of  springs 
and  dashpots  are  often  associated  with  such  integrals  since  they  result  from  the  limit  of  summing 
an  infinite  number  of  Maxwell  and  Kelvin- Voigt  models  (see  [18],  [35]). 

Typical  nonlinear  behaviors  of  the  stress  and  strain  in  rubber  materials  under  finite  ( i.e. ,  non¬ 
infinitesimal)  deformation  include  a  continuous  increase  of  strain  at  decreasing  rates  upon  loading, 
variable  magnitudes  of  the  strain  subject  to  rates  of  loading,  and  different  loading  and  unloading 
paths  due  to  hysteretic  memory  effects.  In  addition,  there  are  other  nonlinear  features  that  are 
particular  to  the  samples  under  study.  These  traits  can  be  modeled  accurately  to  some  degree  using 
theories  for  finite  deformations  alone;  however,  to  fully  describe  nonlinear  viscoelastic  materials, 
it  is  desirable  to  also  consider  internal  chemical  and  physical  interactions  involving  long  chain 
molecules  and  fillers.  We  refer  the  reader  to  [23]  for  the  derivations  of  internal  variable  models  and 
internal  solid  models,  which  are  based  on  the  molecular  point  of  view.  The  previously  mentioned 
Kelvin- Voigt  model,  Maxwell  model  and  Boltzmann  integral  formulation  can  be  modified  to  include 
finite  deformations  by  making  material  coefficients  functions  of  s(t)  and  f,  or  by  defining  new  laws 
between  <j(t )  and  s(t).  There  are  a  number  of  models  that  involve  attempts  to  model  nonlinearity 
in  viscoelastic  material  through  finite  deformation  theories  (see  [18],  [23],  [33] j.  For  more  complete 
discussions  we  refer  readers  to  Chapter  3  of  [35]  and  a  lengthy  list  of  references  found  there. 

Many  of  the  various  models  in  the  literature  have  been  verified  qualitatively  and/or  numerically 
for  certain  samples  under  individual  tests.  However,  due  to  the  complex  dependence  of  rubber 
materials  on  many  physical  parameters,  it  has  been  very  difficult  to  obtain  a  general  formulation 
that  is  reasonably  simple  quantitatively  and  that  captures  viscoelastic  behavior  across  a  wide  range 
of  materials.  For  our  investigations  reported  on  here,  we  have  used  a  Boltzmann  law  with  a  nonlinear 
strain  functional.  Before  turning  to  a  summary  of  our  findings  to  date,  we  briefly  discuss  the  often 
controversial  subject  of  pseudo-phenomenological  models  versus  physics-based  or  internal  models 
in  viscoelasticity. 
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As  we  have  already  noted,  among  the  most  often  used  formulations  to  account  for  hysteresis  is 
the  Boltzmann  superposition  model  (also  called  the  Maxwell  solid  in  history  integral  form)  which  is 
based  on  the  hypothesis  that  the  stress  is  strain  rate  dependent  through  a  convolution  relationship 
[14] 


cr(t)  —  Cs(t)+  f  k(t  -  t)£(t)  dr  .  (16) 

J  —  OO 

For  many  materials,  this  linear  relationship  is  inadequate  to  capture  the  behavior  and  internal 
characteristics  manifested  in  experiments  and  one  must  turn  to  nonlinear  generalizations  of  the 
form 


<r(t)  =  f,  (£(/))  +  [  k(t  —  r)/,.(c:(r))  dr 

J  —  OO 

(17) 

/■*  d 

(/)  =  /,  (.£' ( 1 ) )  +  /  k(t  -  r)— /,  (£•( r))  dr . 

J  —  oo 

(18) 

The  nonlinearities  fe  and  fv  are  often  associated  with  the  material’s  elastic  and  viscoelastic  prop¬ 
erties,  respectively.  Even  if  one  or  both  are  assumed  linear,  they  must  be  identified  along  with 
the  kernel  k  from  experimental  data.  If  this  is  done  with  little  or  no  mechanistic  assumptions 
that  place  constraints  on  the  form  of  /, ,  /,. .  />•  (e.g.,  see  [1] , [2] , [3]  and  references  therein),  such  an 
approach  renders  the  Boltzmann  modeling  a  phenomenological  approach.  If,  on  the  other  hand, 
one  uses  physical  considerations  (such  as  in  [8])  to  constrain  the  choices  of  fe,fv,k  in  the  inverse 
problem  or  parameter  estimation  procedures,  the  modeling  attempt  can  be  thought  of  as  pseudo- 
phenomenological.  In  fact,  it  is  equivalent  to  identifying  the  impulse  response  function  for  certain 
types  of  (possibly  physically-based)  internal  variable  models  (Maxwell  solids  in  differential  form). 
For  example,  if  one  considers  internal  strain  models  as  in  [24],  then  the  basic  assumption  is  that 
one  has  a  finite  number  of  internal  “strain”  variables  £j(t),  j  =  1,  ■  •  ■ ,  N,  along  with  the  strain  s(t), 
and  the  stress  is  given  by 

N 

I  +  W  •  (!9) 

j= i 

The  internal  strains  might  be  given  by  dynamics 

~  =  /?(£«))  »  3  =  !. ' ' '  i  (20) 

or 

^  =  /;(£•(/))£•(/)  ,  j  =  1,  •  ■  • ,  N .  (21) 

If  we  write  the  solutions  of  the  equations  by 

£j(t)=  [  e_(t_T)/Tj/-(e(r))  dr  (22) 

J  —  OO 

or 

£j(0  =  j_x  (23) 

then  (19)  is  completely  equivalent  to  (17)  or  (18)  where 

N 

k(t  -  r)/„(g(r))  =  c^~(t~T)/rjfj(e(T))  (24) 

j= i 
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or 


(25) 


/  N  / 

H*  ~  r)-^fv(£(r))  =  J2  c0~(t~r)/TjJ^  fj«T)) , 

3  = 1 

respectively.  In  this  case,  the  Boltzmann  approach  under  the  assumptions  (24)  or  (25)  is  completely 
equivalent  to  the  internal  variable  approach.  The  only  real  difference  in  these  two  approaches  lies 
in  the  implementation;  i.e.,  one  might  start  with  (17)  or  (18)  and  attempt  to  identify  the  overall 
impulse  response  function  k(t)  without  any  consideration  of  internal  dynamics,  or  one  might  use 
(19)  with  (20)  or  (21)  and  attempt  to  identify  the  individual  decay  constants  Tj  along  with  the 
coefficients  cj  and  fj. 

If,  however,  one  assumes  nonlinear  internal  dynamics  of  the  form 

df  ■ 

:-gf-  +  9j(£j)  =  I,(£j)  ,  j  =  1, ' ' ' ,  N,  (26) 

in  place  of  (20),  then  ( 19)  cannot  be  rewritten  in  Boltzmann  form  and  the  internal  variable  approach 
is  distinct  from  a  general  Boltzmann  approach.  In  this  case,  one  must  solve  a  coupled  system  (the 
overall  dynamical  PDE  in  which  a  appears  plus  the  system  (26)  coupled  with  (19))  as  opposed  to 
an  implicit  internal  dynamics  system  involving  (17)  or  (18)  in  the  PDE’s  of  deformation. 

One  might  also  generalize  the  internal  variable/Boltzmann  linear  models  by  defining  a  gener¬ 
alized  stress  W  =  (a,  a,  ■  ■  ■ ,  cr^”1))  and  generalized  strain  e  =  (e,  £,■■■,  )  with  internal  vector 

dynamics 

—  =  Aa  +  jF(e)  (27) 

so  that 

W(  /  )  ■  •  f  jF(g(r^  dr  . 

J  —OO 

If  we  define  X(t)  =  eAt,  this  becomes 

rt  ni 

/  -  T)J'J(e(T),£(r),---,£(”2)(r))dr,  (28) 

J-oo  j=1 

which  again  is  recognized  as  a  generalized  Boltzmann  formulation.  In  particular,  (18)  is  a  special 
case  of  this  formulation  with  =  1  since 

fv{s(T))  =  K(£{t))£{t)  =  s(  t))  .  (29) 

The  internal  strain  models  of  (19)  and  (20)  can  also  be  put  in  this  more  general  form  by  taking 
£  =  (£i,£2,  •  ■  -  ,£jv)  and 

=  Ae  +  X(s{t)) 

so  that 

e(t)=  f  eA^-T^(£(T))dT 

J  —CO 

and 

ff(i)  =  fe(£{t))  +  c  ■  i(t)  =  /,  (£(/))  +  /  c  •  L  !^-(^.(rj|  (jT  (30) 

for  c  =  (ci,c2,  •  •  ■ ,  cjy ) .  These  considerations  clearly  reveal  that  (19)  and  (20)  are  special  cases  of 
a  generalized  Boltzmann  approach,  and  the  primary  difference  between  a  Boltzmann  and  linear 
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internal  variables  approach  is  mainly  philosophical.  That  is,  in  the  Boltzmann  approach,  one 
attempts  to  identify  the  impulse  response  eAt  as  opposed  to  the  internal  variables  approach  where 
one  is  concerned  with  the  operator  A.  In  both  situations,  one  must  also  identify  the  nonlinearities 
fe  and  fv  (or  T  in  (30)). 

There  is  another  pseudo-phenomenological  approach  to  hysteresis  found  in  the  engineering  liter¬ 
ature  on  viscoelasticity  that  is  related  to  our  remarks  above.  The  so-called  Golla-Hughes-McTavish 
(or  GHM)  method  (see  [20],  [21],  [29]  for  readable  summaries)  can  be  interpreted  as  a  nonphysics- 
based  internal  variable  model  approach  to  viscoelasticity.  In  this  approach,  one  introduces  addi¬ 
tional  coordinates  (i.e. ,  internal  variables)  in  state  space  models  to  account  for  hysteretic  behavior. 
The  general  approach  uses  complex  modulus  or  loss  factor  data  (modulus  or  frequency  data)  to 
fit  rational  polynomials  representing  the  Laplace  transform  of  hysteresis  stress-strain  relationships. 
Specifically,  the  viscoelastic  hysteresis  is  approximated  by  adjoining  a  state  variable  z  with  fre¬ 
quency  domain  representation 

as2  +  f3s 

s 2  +  b.s  +  c 

This  is  equivalent  to  an  internal  dynamics  of  the  form  ave  +  bave  +  cave  =  ae  +  (3i,  or,  in  frequency 
domain,  ove{s2  +  bs  +  c}  =  { as 2  +  / 3s}£ .  This,  of  course,  can  be  recognized  as  a  special  case  of  (27) 
by  defining  W  =  (crve,  <7„e),  e  —  (£,£,£')  and 


a  =  Ko  -f  Ee  (31 ) 

with 

o  n  _  g,o  o  o 

c  b  J  \  0  /3  a 

In  the  usual  engineering  practice,  one  attempts  to  identify  c,  6,  a,  li  using  frequency  domain  data 
and  then  adjoins  this  to  the  time  domain  (e.g.,  finite  element  or  modal)  model  via  (31).  Once  again, 
the  philosophy  of  this  approach  is  to  identify  the  operator  A  (i.e.,  K)  in  (27)  instead  of  identifying 

eAt . 

The  internal  variable  approach  outlined  above  can  be  carried  out  for  particular  structures  and 
materials  by  hypothesizing  physics-based  internal  dynamics.  One  such  example  involves  finite  strain 
internal  stretch  variable  models  for  viscoelastic  rubbers  as  developed  by  [17]  and  [25] .  These  involve 
stick-slip  models  for  the  viscoelastic  response  to  a  large  step-strain  (see  [17]  and  [25]  for  details). 


3.2  Theoretical  Issues 

To  outline  our  efforts  to  date  in  incorporating  hysteresis  in  our  models,  we  return  to  the  example 
involving  simple  extension  in  a  rod  of  cross-sectional  area  Ac,  length  £,  mass  density  p  with  applied 
force  f(t)  at  the  end  x  =  l.  We  assume  a  stress-strain  law  of  the  form 

o{i)  =  £<#(*))  +  /  Y(  t  -  s)-^-gv(£(s),e(s))ds,  (32) 

J  —  oo  (aS 

where  //,  and  are  nonlinear  functions  accounting  for  the  elastic  and  viscoelastic  nonlinear  response 
of  the  rubber  rod.  Our  experimental  quasi-static  data  curves  (see  below)  indicate  that  the  rubber 
follows  different  nonlinear  viscoelastic  constitutive  relationships  in  loading  (e  >  0)  and  unloading 
(e  <  0),  i.e., 


^(£(s),£(s)) 


gvi(s(s))  £(.\l  >  0 
<Jvd(e{  *))  £{s)  <  0, 


(33) 
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where  gvi  and  gv$  are  continuous  nonlinear  functions.  Note  that  //,..,(£(/,))  is  not  necessarily  equal  to 
grcl{£{li ))  at  the  “breakpoints”  fj,  i.e.,  the  points  where  i  changes  sign.  Thus,  we  must  interpret  the 
derivative  in  (32)  in  a  distributional  sense.  Earlier  calculations  and  experiments  suggest  that  the 
rubber  does  not  exhibit  infinite  memory,  but  significantly  depends  only  on  history  of  finite  length 
r.  That  is,  the  memory  kernel  Y  is  such  that  Y(£)  ~  0  for  £  >  r.  Therefore  we  can  approximate 
(32)  by 

°{t)  =  9e{£(t))  +  J  Y(t  -  s)^-gv(e(s),i(s))ds.  (34) 

Let  us  suppose  that  the  t{  are  indexed  so  that  the  rubber  is  unloading  for  //,.  <  /  <  t-k+i,  k  odd. 
Integrating  by  parts  in  (34)  and  assuming  t  /<  <  t  <  tx+ 1,  where  K  is  odd,  we  have 

ft  .  d 

a(t)  =  //,(£(/))  1  j  Y(t-s)—gv(£(s),£(s))ds-Y(r)gvi(e(t-r))  +  Y(0)gvd(e(t)) 

K 

+  ~  tk)(-^)k+1[9m(£{tk))  -  gvd(s(tk))] 

k=  1 

=  //,  (.£(/ >)  f  /  Y(t  -  s)-^gv(£{s),i(s))ds  +  Y(0)gvd(£(t)) 

Jt—r  US 

K 

+  1 );'~ '[//«(£(//.•))  -  gvd(£(tk))\, 

k= 1 


where  we  have  used  Y(r)  &  0.  A  similar  expression  holds  for  loading,  i.e.,  f /<  <  t  <  tx+ 1,  N  even. 
This  is  the  form  of  the  stress-strain  law  that  we  used  in  the  calculations  reported  on  below. 

To  develop  the  corresponding  dynamic  model,  we  must  incorporate  these  jump  terms  into  our 
formulation.  Let  u(t,  x)  denote  the  displacement  at  time  t  of  the  section  of  the  rod  originally  located 
at  X,  0  <  X  <  l.  If  we  assume  that  the  rod  begins  its  motion  at  rest  with  possible  deformation  A{x) 
and  fixed  end  at  x  =  0,  we  obtain  the  model  for  unloading,  i.e.,  tx  <t<  tx+i,  K  odd, 


B^v,  B  (  Bv  Bv  ♦  Bv  B^v 

PA°W  =  Ih  + 

+  Ac  -  tk){-l)k+1[gvi{^{tk))  -  gvd(^{tk))]\  0  <x<l  (35) 

k= 1  J 

(Ac9e(!h)  +  AcY^vd(^)  +  Ac  Jt_r  Y(t  -  S)gv(—{s),  g^(s))ds 
K  f)u  \ 

+4E  -  hX-it+'MgfM)  -  I*-' = /(<}■  (3«) 

In  general  this  model  should  be  written  in  variational  form 

pAlu  -  ±  (am  A)  +  AcY( 0 wa  A)  +  Ac £ Ytt-  dA-ls))d> 

+V  Y,  Y(t  -  lt)(-l}‘+1[Sri(^(l4)  -  SM <  '*))])  =  m  i-  v  (37) 
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for  an  appropriately  chosen  Hilbert  space  V.  This  presumes,  of  course,  that  we  have  sufficient 
smoothness  so  that  evaluation  of  0  at  makes  sense  and  -§^{gvi{  ^ (4)))  G  V* .  Note  that  the 
jump  terms  and  the  related  theoretical  issues  do  not  appear  if  we  suppose  that  gvi  =  gvci  =  gv.  In 
this  case  we  can  take  V  =  if/(0,Q,  H  =  L2( 0,7).  If  we  separate  and  consider  only  the  linear  part, 
i.e.,  if  ge(0  =  E£  +  /?.e(£)  and  gv(£)  =  E£  +  hv(£),  and  we  ignore  he,  hv,  we  obtain 

f)  ft  ♦ 

pAcu  -  ^-(EAC(1  +  y(0))|-  +  EAcj '  Y(t  -  s)^(s)ds)  =  F(t )  .  (38) 

Following  an  approach  we  have  used  in  earlier  treatments  of  Boltzmann  or  time  hysteresis  stress- 
strain  laws  ([1],  [2],  [7],  [8]),  we  may  write  our  model  as  an  abstract  system  in  a  Hilbert  space 
and  develop  a  semigroup  formulation  for  (38)  (see  [7], [8]).  This  in  turn  can  be  used  to  define 
(implicitly)  mild  solutions  for  the  nonlinear  system  in  terms  of  a  nonlinear  variation-of-parameters 
representation. 

For  the  linear  part  of  the  system,  one  can  also  develop  a  variational  formulation  to  obtain  an 
equivalent  well-posedness  framework.  Our  current  efforts  are  devoted  to  developing  the  variational 
formulation  and  corresponding  existence,  uniqueness,  etc.,  for  the  full  nonlinear  system  with  jump 
terms. 

3.3  Experimental  Quasi-static  Results 

In  a  series  of  experiments  at  Lord  Corporation,  we  tested  filled  rubber  rods  in  simple  uniaxial 
tensile  deformations.  Both  ends  of  the  rods  were  manufactured  with  metal  flanges  to  be  secured 
in  a  test  machine,  the  Instron  machine,  which  has  a  movable  load  cell  that  can  apply  tension  to 
one  end  of  the  rod  while  keeping  the  other  end  fixed.  The  software  package  STD4200  was  used  to 
drive  the  Instron  machine  and  to  produce  load- displacement  curves  upon  loading  and  unloading  the 
samples.  Denoting  the  loading  force  by  f(t),  the  displacement  by  A(  ( / ),  the  original  length  of  the 

rod  by  £,  and  the  original  area,  of  the  rod  by  Ac ,  we  calculate  the  engineering  stress  cr(t)  =  and 

the  strain  s(t)  =  X  100%.  Two  types  of  sample  rods  were  used:  ones  filled  with  carbon  black 
(CB)  and  ones  filled  with  silicon  (Sil).  Quasi-static  (u  =  0 ,ii  =  5  in/rnin)  stretching /relaxation 
tensile  cycles  were  performed  on  samples. 

For  each  percent  strain,  the  sample  undergoes  three  cycles  of  loading  and  unloading  to  remove 
possible  small  scale  Mullin’s  effects;  this  results  in  stress  softening  (see  [8]).  We  found  that  the 
first  cycle  had  higher  values  for  the  loads  than  the  second  and  third  loops.  (Mullin’s  effect  is  very 
significant  at  the  initial  puhs  on  the  sample.  In  fact,  the  CB  filled  sample  was  first  pulled  to 
approximately  300%  strain  and  the  Sil  filled  sample  was  puhed  to  150%  strain  before  we  recorded 
data  for  lower  percent  strains.)  We  used  the  third  cycle  for  each  percent  strain  for  our  parameter 
estimation  problems  in  the  results  presented  here. 

We  report  here  on  results  for  a  CB  rod  (£  =  6.13  in,  diameter  =  .54  in  after  removal  of  Mullin’s 
effects)  and  a  Sil  rod  (£  =  7.455  in,  diam  =  .505  in).  We  assumed  no  initial  strain  history  for  the 
rubber  samples,  and  we  used  the  following  stress-strain  relationship  for  our  computations: 

0{t)  =  ge(e(t))  +  j  Y(t  -  s)d9v^e^,e^-ds,  (39) 

Jt-r  US 

where  V  is  the  memory  kernel,  and  <jr  and  gv  represent  the  elastic  response  and  viscoelastic  response 
respectively.  As  before,  r  is  defined  so  that  Y(r)  ~  0. 

This  model  conveniently  provides  a  starting  point  for  modeling  the  hysteresis  curves.  The 
challenge  is  to  identify  plausible  memory  kernels  and  response  functions  ge ,  gv. 

Based  on  previous  experimental  studies  (see  [35]),  we  chose  an  exponential  form  for  the  memory 
kernel  Y.  Such  an  exponential  form  generates  totally  nested  hysteresis  loops  in  the  stress-strain 
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curves,  a  feature  we  also  observed  in  our  experimental  data.  For  C2  and  C3  positive,  we  define 

Y(t)  =  C2e-°3t.  (40) 

The  inclusion  of  both  elastic  and  viscoelastic  response  functions  resulted  from  both  experimental 
observations  and  earlier  work  (e.g.,  see  [23],  [25]).  Elastic  materials  undergoing  finite  deformations 
require  the  use  of  a  nonlinear  constitutive  law,  which  suggests  that  rubber  materials  would  also  have 
a  nonlinear  function  in  the  stress-strain  relationship.  Moreover,  filled  rubber  exhibits  viscoelastic 
and  hysteretic  properties,  which  leads  to  the  use  of  the  viscoelastic  response  function  inside  the 
integral. 

We  tried  a  number  of  linear  and  nonlinear  viscoelastic  response  functions  (see  [35]  for  details). 
The  relative  errors  discussed  in  [35]  suggested  that  a  nonlinear  function  gv(e(t))  is  necessary  for 
both  the  CB  data  and  the  Sil  data  taken  at  sufficiently  large  strains.  Looking  closely  at  the  data, 
it  seemed  that  the  viscoelastic  response  of  the  rubber  during  loading  was  significantly  different 
than  the  response  while  unloading.  As  a  result,  we  include  separate  functions  gvi  and  gvd  in  the 
viscoelastic  response  gv,  where  gvi  represents  the  response  while  e(t )  is  increasing,  and  gv(i  represents 
the  response  when  e(t)  is  decreasing.  This  means  that  gv  is  now  a  function  of  both  e  and  e,  with 
the  form  given  in  (33). 

Additional  curve  fitting  studies  outlined  in  [35]  led  to  our  choice  of  cubic  polynomials  for  the 
viscoelastic  response.  That  is,  gvi(e)  =  J2k=oak£k  and  gvd.(s )  =  J2k=o  where  {«/,-,  bf,}  are 

parameters  to  be  determined  by  estimation  techniques. 

The  form  for  the  elastic  response  ge  was  chosen  in  a  similar  manner.  Preliminary  results 
suggested  that  a  cubic  polynomial  was  also  a  good  choice  for  the  elastic  response.  We  define  ge  by 
ge(s)  =  J2k=o  Ek£k ,  where  Ek  are  parameters  to  be  determined. 

The  parameters  {Ck,  Ek,  a.k,  bk}  were  estimated  by  fitting  stress-strain  data,  i.e.  minimizing 

N 

J2  H**)  -  (41 ) 

k= 1 

for  stress  data  a.  The  single-loop  data  sets  used  for  estimation  include  a  Sil  sample  stretched  to 
100%  strain,  and  a  CB  sample  stretched  to  200%  strain.  These  strain  levels  were  chosen  because 
the  data  exhibited  both  strong  nonlinear  and  hysteretic  characteristics.  The  results  are  presented 
in  Figures  3  and  4.  The  relative  errors  are  at  satisfactory  levels,  and  the  shapes  of  the  loops  are 
closely  matched. 


relative  error  =  2.0% 


strain 


Figure  3:  Sil  rod,  100%  strain:  (left)  stress-strain  with  cubic  ge ,  cubic  gv\  (right)  corresponding 
exponential  kernel. 
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relative  error  =  0.6% 


Figure  4:  CB  rod,  200%  strain:  (left)  stress-strain  with  cubic  ge,  cubic  gv\  (right)  corresponding 
exponential  kernel. 


relative  error  =  0.2% 


strain  time  t 


Figure  5:  CB  rod,  100%  strain:  (left)  stress-strain  with  cubic  ge ,  cubic  gv]  (right)  corresponding 
exponential  kernel. 


After  estimating  the  parameters  with  the  above  data  sets,  we  used  these  parameters  to  predict 
the  stress  for  smaller  strain  loops  of  the  same  rubber  sample.  These  predictions  were  quite  accurate 
for  strain  levels  within  40  percent  of  the  strain  level  used  for  the  original  fitting  ( i.e. ,  strains  at 
180%  and  160%  for  the  200%  strain  data  fit,  and  strains  at  80%  and  60%  for  the  100%  strain 
data  fit),  but  the  accuracy  declined  greatly  as  the  strain  level  decreased.  These  observations  and 
similar  findings  by  [24]  suggest  that  one  parameter  set  may  not  sufficiently  capture  the  stress-strain 
relationship  for  filled  rubber.  Looking  closely  at  the  CB  data,  it  appears  that  both  the  nonlinearity 
and  the  hysteresis  were  quite  different  at  small  strain  levels  versus  large  strain  levels.  Indeed,  these 
observations  led  to  the  conclusion  that  both  a  large  strain  model  and  a  small  strain  model  are 
needed. 

Since  the  parameters  obtained  from  the  CB  200%  data  set  predicted  loops  more  accurately  at 
larger  strains  than  smaller  strains,  we  chose  this  parameter  set  as  our  large  strain  model  for  the 
CB  sample.  We  then  performed  a  separate  fitting  to  the  CB  100%  data  to  obtain  a  small  strain 
model.  These  results  were  similarly  satisfactory  in  terms  of  both  the  relative  error  and  the  shape  of 
the  loop  (see  Figure  5),  and  these  parameters  produced  accurate  predictions  of  the  stress  for  80% 
and  60%  strain  loops. 

In  addition  to  the  accurate  prediction  of  smaller  individual  strain  loops,  another  desirable  feature 
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Silicon:  100%,  80%  and  60%  loops 


Figure  6:  Sil  rod,  response  to  100%,  80%,  60%  strains:  100%  strain  parameter  set,  predicted  nested 
loops  for  80%,  60%  strains;  relative  error  =  3.6%. 


Carbon  black:  200%,  1 80%  and  1 60%  loops 


Figure  7:  CB  rod,  response  to  200%,  180%,  160%  strains:  200%  strain  parameter  set,  predicted 
nested  loops  for  180%,  160%  strains;  relative  error  =  1.2%. 


of  the  stress-strain  relation  is  the  ability  to  predict  nested  loops  consecutively  in  time.  When  the 
parameters  are  used  to  predict  individual  strain  loops,  no  previous  history  is  assumed  and  the 
calculations  are  made  with  one  strain  loop  only.  Predicting  nested  loops,  however,  requires  the 
inclusion  of  all  previous  strain  loops  as  a  part  of  the  history. 

After  obtaining  the  fittings  for  Sil  100%,  CB  200%  and  CB  100%,  the  parameter  sets  were  used 
to  predict  the  inner  loops.  For  the  100%  fittings,  the  inner  loops  correspond  to  80%  and  60%  strain 
levels,  and  for  the  200%  fitting,  the  inner  loops  are  at  180%  and  160%  levels.  Results  are  presented 
in  Figures  6,  7  and  8.  The  predictions  are  quite  satisfactory  with  respect  to  both  relative  error  and 
shape. 
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strain 


Figure  8:  CB  rod,  response  to  100%,  80%,  60%  strains:  100%  strain  parameter  set,  predicted 
nested  loops  for  80%,  60%  strains;  relative  error  =  1.1%. 

4  Concluding  Remarks 

In  the  above  presentation  we  have  outlined  our  progress  to  date  in  the  development  of  nonlinear 
hysteretic  dynamic  models  for  MR  based  elastomers.  Substantial  experimental  validation  for  our 
approach  is  provided  in  the  quasi-static  case.  Significant  efforts  are  still  required  to  achieve  similar 
status  in  the  full  dynamic  case.  Moreover,  to  attain  the  level  of  modeling  required  for  design  and 
feedback  control  of  active  MR  elastomers,  a  major  effort  remains  on  response  of  the  elastomer  in 
the  presence  of  an  external  magnetic  field. 
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